Influence of wetting properties on hydrodynamic boundary conditions at a fluid-solid 

interface 
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It is well known that, at a macroscopic level, the boundary condition for a viscous fluid at a 
solid wall is one of "no-slip". The liquid velocity field vanishes at a fixed solid boundary. In this 
paper, we consider the special case of a liquid that partially wets the solid, i.e. a drop of liquid in 
equilibrium with its vapor on the solid substrate has a finite contact angle. Using extensive Non- 
Equilibrium Molecular Dynamics (NEMD) simulations, we show that when the contact angle is 
large enough, the boundary condition can drastically differ (at a microscopic level) from a "no-slip" 
condition. Slipping lengths exceeding 30 molecular diameters are obtained for a contact angle of 
140 degrees, characteristic of Mercury on glass. On the basis of a Kubo expression for S, we derive 
^ . an expression for the slipping length in terms of equilibrium quantities of the system. The predicted 

behaviour is in very good agreement with the numerical results for the slipping length obtained in 
the NEMD simulations. The existence of large slipping lentgh may have important implications for 
the transport properties in nanoporous media under such "nonwetting" conditions. 



Pacs numbers: 61.20J 68.15 68.45G 



I. INTRODUCTION 



> 

Properties of confined liquids have been the object of constant interest during the last two decades, thanks to the 
considerable development of Surface Force Apparatus (SFA) techniques. While static properties are now rather well 
understood (see e.g. 0), the dynamics of confined systems have been investigated more recently j|. These studies 
have motivated much numerical and theoretical work |^-[Io[| and some progress has been made in giving a simple 
coherent description of the collective dynamics of confined liquids. Both from experimental and theoretical studies 
| has emerged a rather simple description of the dynamics of not too thin liquid films -i.e. films thicker than typically 
10 to 20 atomic sizes. The hydrodynamics of the film can be described by the macroscopic hydrodynamic equations 
with bulk transport coefficients, supplemented by a "no-slip" boundary condition applied in the vicinity (i.e. within 
one molecular layer) of the solid wall. Hence, in spite of the fact that the wall induces a structuration of the fluid 
into layers that can extend 5 to 6 molecular diameters from the wall, the hydrodynamic properties of the interface 
are quite simple. 

It turns out, however, that all experimental and numerical studies of confined fluids that have been carried out with 
fluid/substrate combinations that correspond to a total wetting situation. By this we mean that 7ls + Jlv < Jsv, 
where 7 is a surface energy and the indexes L,V and S refere to the liquid, its vapor and the substrate, respectively 
p7[ . In this letter, we investigate the structure and the hydrodynamic properties of a fluid film that is forced to 
penetrate a narrow pore in a situation of partial wetting, i.e. when jls + Jlv > Isv ■ This corresponds to the case 



where a drop of the liquid resting on the same substrate, at equilibrium with its vapor, has a finite contact angle 
which can be deduced from Young's law ^lv cos 6 = (jsv — Ils) Mi- 



ll. MODEL AND RESULTS 



We first describe our model for the fluid and the substrate, and some details of the simulation procedure. All 
interactions are of the Lennard- Jones type, 

%(r)=4e((^) 12 -c ii (^) 6 ) (1) 
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with identical interaction energies and molecular diameters a. The surface energies will be adjusted by tuning the 
coefficients Cy . In all the simulations that are presented in this letter, the solid substrate is described by atoms fixed 
on a FCC lattice, with a reduced density psa 3 = 0.9. As the atoms are fixed, the coefficient ess is in fact irrelevant. 
The interactions between fluid atoms are characterized by cff = 1.2, meaning that the fluid under study is more 
cohesive than the usual Lennard- Jones fluid. The fluid-substrate interaction coefficient cfs will be varied between 
0.5 and 1. All the simulations will be carried out at a constant reduced temperature, ksT/e = 1. 

Finally, we mention that the configuration under study will be that of an fluid slab confined between two parallel 
solid walls. Typically, a configuration contains 10000 atoms, with a distance between solid walls h = 18a and lateral 
cell dimensions L x = L y = 20a. Periodic boundary conditions are applied in the x and y directions, i.e. parallel to 
the walls. For each wall, three layers of FCC solid (in the 100 orientation) will be modelled using point atoms, a 
continuous attraction between the fluid and the wall in the direction perpendicular to the walls being added in order 
to model the influence of the deeper solid layers. 

All simulations were carried out at constant temperature (fceT/e = 1) by coupling the fluid atoms to a Hoover's 



thermostat 15 1. In flow experiments, only the velocity component in the direction orthogonal to the flow was 
thermostatted. 

Before we discuss in detail the structure of a film, we can roughly estimate the influence of the interaction parameters 
on the wetting properties of the fluid. Following the standard Laplace estimate of surface energies [Q, we have 
7ij = —PiPj J™ ruij(r)dr. Using Young's law, we obtain for the contact angle cos# = — 1 + 2psCFs/ Pfcff- 

From this expression a variation of cfs between 0.5 and 0.9 would be expected to induce a variation of 9 between 
100 degrees and 50 degrees. A more accurate determination of the surface tensions was carried out using the method 
of Nijmeijer et al. The surface tensions are defined in terms of an integral over components of the pressure tensors 
which can be computed in a simulation. We refer to ref. jl6| for more details. The results are listed in table I. By 
tuning the solid-fluid interaction strength cfs from cfs = 1-0 to cfs = 0.5, we found the contact angle deduced from 
Young's law varies from 9 ~ 90 degrees to 9 ~ 140 degrees. In figure |], a typical configuration of a liquid droplet (in 
coexistence with its vapor) on the solid substrate is shown for cps = 0.5, corresponding to 9 ~ 140 degrees. In the 
following we shall describe such large contact angles (i.e. larger than 90 degrees) as corresponding to a "nonwetting" 
situation. 

In order to force the fluid into a narrow liquid pore under such partial wetting conditions, an external pressure has 
to be applied. A simple thermodynamic argument shows that for a parallel slit of width h, the minimal pressure is 
■Po = 2(lLS-lSv)/h. For the fluid with c FS = 0.5, we find, for h = 18a, P a = 0.079e/cr 3 , while P = 0.018e/cr 3 when 
cps — 0.9. If we use a — 5A, e = 0.05eV, then Pq ~MPa for h = 9nm in the "nonwetting" case, 9 — 140 degrees. 

Figure |^ shows the density profiles of the nonwetting fluid inside the pore for pressures corresponding to 2.8-Po and 
16.4Po- The pressure is changed at constant pore width by changing the number of particles. It is seen in this figures 
that the highest pressure structure strongly resembles what would be obtained for the usual case of a wetting fluid, 
with a strong layering at the wall. The structure at the lower pressure is markedly different, with both a layering 
parallel to the wall and a density depletion near the wall. 

We now turn to the study of the dynamical properties of the confined fluid layer. Two types of numerical experi- 
ments, corresponding to Couette and Poiseuille planar flows were carried out. In the Couette flow experiments, the 
upper wall is moved with a velocity U (typically U = 0.5 in reduced Lennard- Jones units). In the Poiseuille flow 
experiment, an external force in the x direction is applied to the fluid particles. In figures Q and ^ we compare the 
resulting velocity profiles to those that would be expected for a " no slip" boundary condition applied at one molecular 
layer from the solid wall. Obviously the velocity profiles for the nonwetting fluid imply a large amount of slip at the 
solid boundary. As usual, this slippage effect can be quantified by introducing a "partial slip" boundary condition for 
the tangential velocity vt at the solid liquid boundary: 

dv t{ 1 , ,„s 

^■U =aw =-^U„. (2) 

This boundary condition depends on two parameters, the wall location z w and the slipping length 5. By studying 
simultaneously Couette and Poiseuille flow for the same fluid film, both parameters can be determined if they are used 
as fit parameters for the velocity profiles obtained in the simulation. The results of such an adjustment are shown in 
figures U and [|. It turns out that, as was the case in earlier studies ||, the hydrodynamic position z w of the walls 
is located inside the fluid, typically one atomic distance from the outer layer of solid atoms. Much more interesting 
is the variation of the slipping length 5, which in earlier work was always found to be very small. In figure the 
variation of S as a function of the pressure is shown for several values of the interaction parameters. The pressure 
are normalized by the capillary pressure Pq defined above as the minimal pressure that must be applied to the fluid 
in order to enter the pore. For an interaction parameter cfs = 0.9, corresponding to a contact angle 9 — lOOdegrees 
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, the usual behaviour (i.e. a small 5) is obtained. For an interaction parameter cps — 0.5 which corresponds to a 
contact angle 9 = 150 degrees p8| , slipping lengths larger than 15 molecular diameters can be obtained at the lowest 
pressures; even at relatively high pressures (10Po)j the slipping length remains appreciably larger than the molecular 
size a. 



III. THEORY 



In order to understand the relation between the hydrodynamic boundary condition and the wetting properties of 
the fluid on the substrate, one needs to estimate the dependence of the slipping length S on the microscopic parameters 
of the system, such as "roughness", temperature, cfs, etc... 

From the "kinetic" point of view, this is obviously a hard task to complete since the slipping length accounts for 
the 'parallel' transfer of momentum between the fluid and the substrate. Explicit calculations can only be done in 
some model systems |l9| . Now in the general case of a dense fluid, no explicit formula for 5 in terms of microscopic 
quantities is available in the litterature. 

In the following, we derive an approximate expression for the slipping length, which will allow us to discuss quali- 
tatively the relationship between wetting and boundary conditions. 

Our starting point is a Green-Kubo expression for the slipping length || : 

A = ? = TTT^r ^ (F x (t).F x (0)) dt (3) 
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where rj is the shear viscosity of the fluid, A — L x L y is the lateral solid surface, and F x is the x-component of the 
total force due to the wall acting on the fluid, at equilibrium (x is a component parallel to the wall). The quantity 
A = rj/8 can be interpreted as the friction coefficient of the fluid/wall interface, relating the force along x due to the 
wall, to the fluid velocity slip at the wall 

(F x ) = -A\ V sHp (4) 
By introducing the density-density correlation function, eq. (0) can be rewritten 
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l — £ dtj dn J dr 2 F x (n) F x (r 2 ) (p(n;t) p(r 2 ;0))> 



(5) 



The force field F x (r\) derives from the fluid-substrate potential energy V(xyz). The latter has been computed by 
Steele for a periodic substrate interacting with the fluid through Lennard- Jones interactions pofl . In the case of a 
(100) face, the main contribution can be written in terms of the shortest reciprocal lattice vectors according to 

V(xyz) = V (z) + Vi(z) {cos (q^x) + cos (q\\y)} (6) 

where qu = 2n/a s and a s is the lattice spacing in the fee solid. The altitude z is the distance to the first layer of the 
atoms of the solid. The functions Vq(z) and V\{z) are given by Steele in ref ( pQ] ). In our case however, the attractive 
parts of Vq(z) and Vi(z) are multiplied by the tuning factor cps- The force along x is the derivative with respect to 
x of the interaction potential V(xyz) : 

F x (x,y,z) = q\\ Vi(z) sin(g||x) (7) 

When inserted into eq. (||), one obtains 
r\ 1 



. j; T j dtj dxidyidzi j dx 2 dy 2 dz 2 <jjjVi(zi)Vi (z 2 ) sin (t?||a;i) sin (q\\x 2 ) (p(ri;t) p(r 2 ; 0))) (8) 
We now introduce the Fourier transform of the density in the plane parallel to the substrate 

Pk(z)(t) = J dxp( X ,z;t) e i k x (9) 
with x = (x, y) and the vector k is parallel to the substrate. This allows to rewrite the previous equation (pj) as 
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dzi / dz 2 Vl(zi) Vx(z 2 ) / dt»((p B|1 («i)(t)^- au («2)(0)» (10) 







with gii in the x direction and -ft stands for the real part. Note that in deriving eq. (|Io|), the homogeneity of the 
system in the direction parallel to the substrate was taken into account. 

Due to the presence of the confining solid, we now assume that the main contribution of C(q\\, z\, z 2 ;t) to the 
time-integral comes from the dynamics in the plane (x, y) parallel to the substrate. Moreover, these dynamics are 
probed at the first reciprocal lattice vector qn. Since q\\ is close to the position of the first peak in the structure factor, 
it is reasonable to assume a diffusive relaxation of (p q (zx)(t)p- q , ] (z 2 )(0)) |pi|)p2[ , yielding 

(Pq {l (2i)(*)P-«n ( z 2)(0)) = exp(-q\ D n t) {p n (zi)p- n (z 2 )) (11) 

where D qil is a collective diffusion coefficient and (p qil (zi)/3_ 9|| (^2)) is the static correlation function. 
The time integration can now be performed to obtain 
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dzi I dz 2 Vi(zi) Vi(z 2 ) (p q {zi)p- q (z 2 )) (12) 



so that the slipping length is expressed in terms of static properties of the inhomogeneous system only. A further 
simplifications can be done by assuming that, due to the stratification near the substrate, the main contribution in 
the integrals in eq. ( |l2|) arises from the z\ ~ z 2 terms, so that {p qn (zi)p^ qn (z 2 )) ~ A(p(zi)) 5{z\ — z 2 ) S(q\\ \z\). 
The quantity 5(q|||zi) is the z dependent structure factor, in the plane z = z\ (parallel to the solid), defined as 

The factor A{p(zi)} (A being the lateral surface) normalizes the average by the number of fluid particles in the layer 
at the altitude z\. If no locking of the fluid occurs near the substrate, one may approximate ^(gy 1 21 ) by its value at 
the first layer, S(q^\z 1 ) ~ Si(qn). Equation ( [l2| ) thus reduces to 

dz x p{z x ) V^z,) 2 (14) 
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In order to have a practical estimate to compare with, this formula can be further approximated. The integral term 
in eq. ( |l4| ) may be approximated by assuming that it is dominated by the behaviour around the first layer located at 
z c ~ a, so that 

f r + co 

dz p(z) V x (zf ~ p c / dz Vi{z) 2 (15) 

J a 

with p c the density at the first layer, denoted in the latter as the "contact" density. Moreover, one expects the "long 
range" attractive part of Vi(z) to contribute mainly to the integral. Since the latter can be written cfs Vi tt (z), with 
Vi tt (z) independent of cfs, one gets 



a Si{q\\) c% s p c a 3 



« 5 s J" I .a (16) 



where D* — D qil /Do, and Dq — k B T /2>nr]a is the Stokes-Einstein estimate for the bulk self diffusion coefficient. 

All the quantities involved in eq. ( |l6| ) can be computed in equilibrium Molecular Dynamics simulations. The 
density at contact p c can be measured from density profiles, such as in fig. ^. On the other hand, Dq., and 
Si(q\\) can be computed from the correlations of density fluctuations in the first layer. In practice, we introduce 
the function Si(q\\,t) = N^ 1 (p q {t)p q ^ (0)), where N\ is the average number of particles in the first layer and 
Pq\\tf) = J2k=i N 1 ex Pi c l\\ Xk (t) is restricted to atoms in the first liquid layer. The value of Si(q\\,t) at time t = 
yields S±(q\\), while D q][ is obtained in terms of the inverse relaxation time of Si(q\\,t), according to eq. ( |TT|) . Let us 
note at this point that the assumption of an exponential decay of Si(q\\,t) is indeed verified in the simulations, which 
allows us to clearly define Z? 9|| . 

In fig. ||, the ratio 8/5*, with S*/a = D* / Si(q\\)c FS is plotted as a function of the inverse density at contact, 
l/p c a 3 . In these variables, the theoretical estimate, eq. (pi), predicts a linear dependence of 8/5* as a function of 
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l/p c a 3 . As shown in fig. ||, a linear behaviour 5/8* = a(l/p c cr a — 1/ p s hiftO~ 3 ) is indeed observed, in agreement with 
the prediction. A least-square fit of the datas in this plot gives a = 3.04 and l/p s hift& 3 = 0.47. The presence of a 
shift in the density, l/p s hift& 3 , can be interpreted to account for the higher-order correction in the density at contact 
which have been neglected in deriving eq. (|l7J) (in particular in the rough approximation assumed in eq. (|l5|)). In 
the interesting limit where the contact density p c is small and 5 is large, this shift does not contribute anymore. 
In hg. ||, the full theoretical result for 6 
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a Si(q\\) c FS p c a 3 \ p B hift 

is plotted as function of P/Po against the measured (out-of-equilibrium) results for 6. The good agreement obtained 
in these variables for all different pressures and interaction strength cfs emphasize the robustness of the previous 
theoretical estimate. Obviously, this expression breaks down for very large contact density rho c > rho s hift = 2.1cr 3 , 
where d is expected to vanish anyway. 

This result calls for several comments. First, eq. (|l7j ) shows that, at for given fluid-substrate interaction cfs, 
S decreases with the density and structuration of the fluid in the first layer. The slipping length is thus expected 
to be quite small in a dense fluid at high pressures, as usually observed and measured experimentally [8|JT^]. More 
specifically, eq. ( |l7j ) predicts a strong dependence of S on the value of the structure factor in the first layer, taken at 
the shortest reciprocal lattice vector qu. This result is in qualitative agreement with previous simulation results 
Now if the fluid-substrate interaction cps is decreased at a given contact density of the fluid p c (eg, by increasing 
simultaneously the pressure), the previous result predicts a strong increase of the slipping length. This explains 
why substantial slip may be obtained, even if a strong structuration does exist in the fluid, a fact which is a priori 
counter- intuitive. 

Finally, let us come back to the problem of the influence of the wetting properties on the slipping length 5. As 
noted above, eq. ( $7\ ) predicts that 6 is a decreasing function of the interaction strength cfs- Now, as emphasized 
for example by the Laplace estimate of the contact angle, cos# = —1 + 2p$CFS I ' PfCff, the contact angle may be 
interpreted as a "measure" of the fluid-substrate interaction strength cfs- I n particular one expects the fluid to 
approach a non- wetting situation (cos# — > —1), when cfs decreases to zero. The previous equation, eq. (|l7|), thus 
predicts a strong increase of the slipping length S when cos 8 — > — 1. In other words in the idealized situation of a non 
wetting fluid, 9 = it, a perfect slip may be expected for the boundary condition of the fluid near the surface. The 
correct trend is observed in our simulation results. This result is in agreement with several experimental observations 
pjyjfl, reporting very large slipping lengths for nonwetting fluids. 



IV. CONCLUSIONS 



Obviously the existence of such a large slippage effect should manifest itself in the dynamical properties of a liquid 
confined in a nanoporous medium. If one considers a single cylindrical capillary, a straightforward calculation in the 
Poiseuille geometry shows that the existence of slip on the boundaries increases the flow rate in the tube as compared 
to the "usual" no-slip case by a factor 1 + 86 /h (with h the pore diameter and 5 the slip length). Thus, in a porous 
medium, the effective permeability K e ff, which relates according to Darcy's law the flow rate to the pressure drop 
p3| , is expected to increase by the same factor : 

K ef f=K (l + 8~) (18) 

where Kq is the "standard" permeability, obtained within the no slip assumption (i.e., when 6 is zero). In a wetting 
situation, 5 is obtained to be very small and K e ff ~ Kq. However, in a nonwetting situation (6 ~ 140 degrees), the 
slipping length 6 may largely exceed the nanometric pore sizes h, so that the effective permeability K e ff is expected 
to be much larger than Kq (say, more that one order of magnitude in view of the prefactors) . 

It can also be expected that the microscopic dynamics of the molecules could be rather different in a " nonwetting" 
medium, compared to what it is in the bulk or in a medium with strong solid/liquid affinity. In fact, recent studies 
point towards the importance of the surface treatment for the reorientation dynamics of small molecules in nanopores 
. Correlating the wetting properties with such microscopic studies seems to be a promising area for future research. 
This work was supported by the Pole Scientifique de Modelisation Numerique at ENS-Lyon, the CDCSP at the 
University of Lyon, the DGA and the French Ministry of Education under contract 98/1776. We would like to thank 
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TABLE I. Dependence of the surface tensions (in units of e) and contact angle 9 (in degrees), on the tuning parameter cps- 
The liquid-vapour surface tension was determined independently to be 'Ylv = 0.94 e. 



CFS 


7s v — JLS 


cos(6>) 


9 


0.5 


-0.71 


-0.74 


137 


0.6 


-0.65 


-0.68 


133 


0.7 


-0.50 


-0.52 


121 


0.8 


-0.35 


-0.36 


111 


0.9 


-0.16 


-0.17 


99 



G 



o 




FIG. 1. A typical configuration of a liquid droplet (1000 atoms) on a solid substrate in a "nonwetting" case (cfs = 0.5), 
in equilibrium with its vapor. This configuration was prepared by starting from an homogeneous fluid slab of thickness 18<r, 
confined between two walls. The droplet is formed by simultaneously increasing moving the upper wall by 20a in the z direction 
and removing the fluid atoms that lie near the box boundaries. The bounding box has a size of 20a. 
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FIG. 2. Density profiles of the "nonwetting" fluid (cfs = 0.5) confined between two solid walls separated by 20<r. The 
positions of the first layer of solid atoms have been indicated by vertical dashlines. Full line: P/Po = 2.8; dashed line: 
P/Po = 16.4. The latter curve has been shifted upwards for clarity. 
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FIG. 3. Velocity profile (in reduced units) of the "nonwetting" fluid (cfs = 0.5) in a Couette geometry. The reduced pressure 
is P/Po — 7.3. The solid line is the simulation result, the dashed line is a linear fit of the numerical results, and the dashed-dot 
line is the velocity profile predicted by the no-slip bounbary condition. The velocity of the upper wall is U = 0.5. 




FIG. 4. Same as in figure 3, but in a Poiseuille geometry. The external force applied on each fluid particle is f ex t = 0.02e/a 
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FIG. 5. Variation of the slipping length 8 (in units of a) as a function of the reduced pressure P/Po, for several values of 
the interaction parameters cfs- From top to bottom, the datas correspond to cfs = 0.5, cfs = 0.6, cfs ~ 0.7, cfs = 0.8, 
cfs = 0.9. Solid lines are the theoretical prediction, eq. ( |l7| ) (see text for details). 




FIG. 6. Normalized slipping length 8/8* , with 8* — cD*. /Si(q\\)c F s, as a function of the inverse contact density l/p c cr 3 . In 
this plot, a linear dependence is expected according to the theoretical prediction, eq. (^). The dashed line is a least-square fit 
of the numerical datas, with slope a = 3.04 and shift in inverse density \j p s hift<J A = 0.47. 
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